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ABSTRACT 

We study axisymmetric models of layered protoplanetary discs taking radiative trans¬ 
fer effects into account, and allowing for a residual viscosity in the dead zone. We also 
explore the effect of different viscosity prescriptions. In addition to the ring instability 
reported in the first paper of the series we find an oscillatory instability of the dead 
zone, accompanied by variations of the accretion rate onto the central star. We pro¬ 
vide a simplified analytical description explaining the mechanism of the oscillations. 
Finally, we find that the residual viscosity enables stationary accretion in large regions 
of layered discs. Based on results obtained with the help of a simple 1-D hydrocode 
we identify these regions, and discuss conditions in which layered discs can give rise 
to FU Orionis phenomena. 
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1 INTRODUCTION 

The transport of angular momentum in protoplanetary discs 
is most probabl y driven by the magnetoro tational instabil¬ 
ity, described bv iBalbus fc Hawlevl (Il991al ral. and hereafter 
referred to as MRI. Since the protoplanetary discs are ion¬ 
ized rather weakly, in some regions of the disc the ionization 
degree ^ can decrease below a critical level, ^c, at which the 
gas decouples from the magnetic field and the MRI decays. 
Such region is referred to as a dead zone, contrary to the so 
called active zone where the MRI operates. 

The most important processes responsible for the ion¬ 
ization of gas in protoplanetary discs are 1 ) particle colli¬ 
sions due to thermal motions, and 2 ) irradiation by cosmic 
rays and high-energy photons originating from the central 
star. The collisional ionization dominates in the innermost 
part of the disc, where the temperature is high enough to 
maintain ^ at all distances from the disc plane. Further 
away from the star is exceeded only in two surface layers 
(on both sides of the disc), each ionized by cosmic rays and 
stellar X-ray quanta, and having approximately constant 
column density Ea. Both processes result in the following 
distribution of disc activity: at r < ri (where the thermal 
ionization degree exceeds ^c) the whole disc is active; at 
ri < r < r 2 (where E > Ea) a dead zone is sandwiched be- 
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tween two active layers; and at r > r 2 (where E < Ea) the 
w hole disc become s active again. This model was proposed 
by|G ammid Jl99(f> to explain the FU Ori outbursts. 

FU Ori stars are young stellar objects whose optical 
brightness occasionally increases by 3—4 mag on a time-scale 
of 1 — 10 yr, and remains at the increased level for 10 — 100 yr. 
It is generally believed that these events are related to the 
variations of the accretion rate in circumstellar discs. To 
explain the observed increase in luminosities the accretion 
rate has to change from ~ 10 ~^ to ~ 10 ~^ M 0 yr“^, implying 
that as much as IO^^Mq must be accreted onto the central 
star during the outburst iHartmann fc KenvorJll99d) . 

The FU Ori phenomenon may be caused by a 
thermal instabili t y of t he protoplanetary disc (see e.g. 
iLodato fc Clark3 JioO^), iBell fc LirJ (Il994fl and references 
therein), analogous to the one operating in accretion 
discs around primary components of cataclysmic binaries 
iMever fc Mever-Hofmeisted Il98lh . However, the thermal 
instability model has difficulties with explaining long du¬ 
ration of the outbursts, since the high accretion rate dur¬ 
ing the outburst quickly removes most of the material from 
the inner part of the disc where the instability occurs 
iHartmann fc KenvorJlT99(il . A transition of the very inner 
disc to the outburst sta t e due to the thermal instability was 
studied bv iKlev fc LinI Jl999ll with the help of the similar 
two-dimensional code as used in this work. 

lOarnmi^ il99(dl suggested that this problem may be 
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solved by the layered model, in which the dead zone serves 
as a mass reservoir for the outburst. In the layered disc the 
accretion rate can be an increasing function of r. When this 
is the case, an annulus centred at ro receives more mass 
per unit time from the disc exterior to it (r > ro than it 
loses to the disc interior to it (r < ro), causing the matter 
to accumulate in the dead zone. If this process lasts long 
enough, the total surface density of the disc may become 
so large that once the accretion was somehow triggered in 
the dead zone it would release enough heat to restore the 
coupling between the matter and the magnetic field. Such 
an ’’ignition” of the accumulated matter would then start a 
self-sustaining outflow of mas from the dead zone into the 
inn er disc, substantially in creasing the accretion rate. 

lArmitage et al.l assume that the triggering fac¬ 

tor is the gravitational instability. Their model shows re¬ 
peating periods of vigorous accretion separated by quies¬ 
cent intervals lasting ~ 10® yr, in qualitative agreement 
with the observed properties of FU Ori objects. However, 
our radiation-hydrodynamic simu lations of layered discs 
iWiinsch. Klahr fc R,6zvczkall2nn^ hereafter Paper I) indi¬ 
cate, that it may be rather difficult for the dead zone to store 
~ 1O“^M0 of the disc matter necessary to feed an FU Ori 
outburst. The reason is the ring instability described in Pa¬ 
per I, which, long before the required mass is reached, causes 
the dead zone to split into well-defined rings that become un¬ 
stable to the instability discovered bv IPapaloizou fc Pringlel 

JTgil. 

In the present paper we extend the simulations reported 
in Paper I onto a broader range of disc parameters, with the 
principal aim to hud in which conditions the dead zone could 
become massive enough to account for FU Ori outbursts. 
For reasons explained in §2, along with models whose dead 
zones are entirely free of viscosity we explore models with 
some residual viscosity acting in the region which formally 
should be declared dead. We also explore the effect of dif¬ 
ferent viscosity prescriptions. 

The paper is organized as follows. In §2 we briefly de¬ 
scribe the numerical methods and input physics. The results 
of simulations are reported in §3. In §4 we give a simple an¬ 
alytical description of processes at the inner boundary of 
the dead zone, and derive conditions for ignition. In §5 we 
summarize our results and discuss the limits of the mass 
accumulation process. 


2 NUMERICAL METHODS AND INPUT 

PHYSICS 

As in Paper I, the simulations are performed in spherical 
coordinates wit h the help of the radiation-h ydrodynamics 
code TRAMP (iKlahr. Henning fc Klevlll99flll . The models 
are axially symmetric, but we do not impose an additional 
mirror symmetry with respect to the mid-plane of the disc, 
i.e. we allow for the flow through the mid-plane. Further 
technical details, like boundary conditions or initialization 
procedure, can be found in Paper I. 

The viscosity is treated a ccording to the a-prescription 
of IShakura fc SunavevI (Il973ll . In Paper I, following the ar¬ 
gument that the turbulent viscosity coefficient should be 
proportional to the product of sound speed and the largest 
scale available for the turbulent eddies, we used the formula 


V = acB,iHi, , ( 1 ) 

where is the half-thickness of the active layer, Cs,i is the 
sound-speed at the bottom of the active layer and q is a di¬ 
mensionless parameter discussed below. In the present work 
this formula is employed in models 1 — 3. The remaining 
models are obtained based on an alternative formula which 
allows for variations of i/ with the distance from the mid¬ 
plane 



where Cs and Q are local values of sound speed and angular 
velocity. Obviously, in a vertically averaged approximation 
assuming the zero-thickness dead zone both prescriptions 
are equivalent; however, in two-dimensional models with the 
dead zone, we found that the details of the evolution are 
sensitive to the particular form of the viscosity formula (see 

§ 3 ). 

We assume that at T > Tuni = 1000 K the ionization 
degree is high enough for the magnetorotational instability 
to operate and provide the source of viscosity, whereas at 
T < Tiim the MRI is quenched. This assumption is based 
on the fact that at T ^ 1000 K ^ increases b y hve orders 
of m agnitude due to ionization of potassium dUmebavashil 
ll 98,3^ . Furt her, based on the cosmic ray s topping depth cal¬ 
culated by lUmebavashi fc Naka^ (Il98Ji . we assume that 
on both sides of the disc a surface layer of column density 
Ea = 100 gcm“^ is made active due to ionization by ener¬ 
getic particles of cosmic origin. 

In principle, ionization due to X-rays from the central 
star should also be taken into account. However, the corre¬ 
sponding active surfa ce layer is probably not thicker than 
a few tens of gcm“^ llOlassgold. Naiita fc Igealll997^ . The 
X-rays become important when the disc i s screened from 
cosm i c rays by magnetized stellar winds llGlassgold et alj 
ll 997t iFromang et all l2^^2^ . In such cases the ionization 
structure of the disc strongly depends on various model 
parameters like properties of the X-ray source, critical 
mag netic Reynolds num ber or abundance of heavy met¬ 
als llFromang et a, Dl», and it is not possible to define 
one critical Ea for all possible circumstances. Unfortunately, 
such effects are too complicated to take them into account 
in multidimensional simulations, and an active surface layer 
of constant column density must serve as the necessary ap¬ 
proximation. In effect, to define the boundary of the dead 
zone we emplo y the simple mo del of disc ionization structure 
introduced bv iGammia Jl99(f) . 

In Paper I a was set to 0 within the dead zone. However, 
it has been recently argued that the dead zone is not entirely 
dead in the sense that a residual viscosity can be main¬ 
tained there by a purely hydrodynamic turbulence, excited 
by the non-axisymmetri c waves propagating fro m the MRI- 
turbulent active layers (iFleming fc Stoii3l2003ll . Following 
those arguments, in some models we allow for odz A 0. 
However, in all models the dead zone is much less viscous 
than the active zone, i.e. aoz Oa- 

In the previous work we neglected the transport of heat 
by micro-turbulence, i.e. the heat was transported across the 
disc by radiation only. The same approximation is adopted 
here. W hile it would b e poorly justified in 1-D models 
(see e.g. iGa.nniz^ l2nnih . it is entirely reasonable in 2-D 
discs, where the most important (i.e. the largest) turbu- 
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Model 

jy-type 


odz 

grid 

1 

aCuHa 

10-2 

0 

128 X 32 

2 

aCsHa 

10-2 

10--^ 

128 X 32 

3 

aCsHa 

10-2 

10-3 

128 X 32 

4 

aCg /n 

10-2 

0 

128 X 32 

5 

QCg /n 

10-2 

10-4 

128 X 32 

6 

ac^/n 

10-2 

10-3 

256 X 64 

7 

oCg/U 

5 x 10-3 

5 X 10-4 

128 X 32 

8 

ac'^/n 

2 X 10-2 

2 X 10-3 

128 X 32 


Table 1. Parameters of the disc models. 


lent scales are automatically taken into account, as they are 
resolved on the grid. Based on the similarity of the tur¬ 
bulence maintained by th e magnetorotational instability to 
Kolmogoroff turb ulence (iHawlev. Gammie fc Balbujll99^ 
IStone et alJll99^ . we estimated that the subgrid micro- 
turbulent flux is not significant compared to the radiative 
one. 


3 RESULTS 

Varying a^, odz and the viscosity prescription, we obtained 
a set of disc models listed in Table Q Models 1-5, 7 and 8 
were followed for several thousand orbits of the outer edge of 
the disc on a grid of 128 x 32 points in r and 9, respectively. 
Model 6 , which is compared to the analytical description of 
the disc obtained in is computed on a grid with double 
resolution (256 x 64). The size of the computational domain 
was set so that it included a small part of the inner a-disc 
and a part of the dead zone where the ignitions took place. 
In radial direction, it was r = (0.05, 0.35) AU for models 
1-6, r = (0.04, 0.3) AU for model 7 and r = (0.1, 0.7) AU for 
model 8 . The vertical extent of the computational domain 
was 6 = (—5°, 5°) for all models. 

3.1 Model 1. {v = aCsHa, eta. = 10“^, aoz = 0) 

This model is the same as the one described in Paper I, 
but the current resolution is two times lower. The main 
purpose of this calculation was to check if we recover all 
essential features and time-scales observed and measured 
in Paper I, which is indeed the case. The model develops 
a strong instability, which quickly leads to the formation 
of well-separated rings, most of which remain stable for at 
least a few hundred years. However, they would probably 
de cay in 3D due to the hydrody namic instability discoyered 
by IPaualoizou fc Pringla ^ll985^ as noted in Paper I. 

Paper I contains a detailed discussion of the ring in¬ 
stability mechanism. Here we only mention that the key ele¬ 
ment of that mechanism is the coupling between the changes 
of the dead zone thickness and the mass accumulation rate. 
It turns out that the accretion rate drops at the inner edge 
of the ring, which works as a bottle-neck in the accretion 
flow along the active surface layer. As a consequence, mat¬ 
ter accumulates preferably in the ring. 

Since only a small amount of mass flows from the ring 
inwards, the growth of a given ring is substantially slowed 
down when additional rings are formed further away from 
the centre of the disc. If several rings are simultaneously 
present in the grid, the accretion rate (always measured at 


the inner boundary of the computational domain) drops to 
a value as low as ~ 3 x 10“® M 0 yr“^. As the ring accretes 
mass and becomes denser, its central temperature increases. 
Once it exceeds Tiim, the ring ignites and a ’’mini-outburst” 
ensues, after which the ring is rebuilt at the same position. 
The accretion rate onto the star increases up to to approx¬ 
imately 5 X 10“® M0yr“^ during an ignition of a single in¬ 
nermost ring (at t = 80 and 220 yr). Moreover, two cases of 
’’induced ignition” are observed at t = 490 and 680 yr. In 
such a case the heat produced by the ignition and accretion 
of a particular ring is transported by the radiation to the 
neighbour rings and makes them active. The accretion rate 
onto the star increases up to 1.4 and 1.0 x 10“® M 0 yr“^. 


3.2 Model 2. {v = aCsHa, eta = 10 odz = 10 "') 

This model evolves very similarly to Model 1. The non-zero 
viscosity in the dead zone causes the temperature inside the 
rings to increase faster, and the induced ignition occurs ear¬ 
lier than in Model 1. The accretion rate in both quiet and 
ignited states are almost the same as in Model 1. 


3.3 Model 3. {u = aCsHa, eta = 10 odz = 10 

In this model we also observe the tendency toward ring for¬ 
mation. However, the viscosity in the dead zone is now suf¬ 
ficiently high to prevent it from splitting. The behaviour 
characteristic of the rings in Models 1 and 2 is now observed 
in the whole low-viscosity area: as the surface density of the 
disc grows, the temperature at the mid-plane grows as well, 
and ignition occurs once it exceeds Tum. The accretion rate 
is increased to ~ 1.3 x 10“® M 0 yr“^, the accumulated mass 
is quickly removed from the ignited region, and the layered 
structure is re-established. Such mini-outbursts repeat on 
a time-scale of 100 — 200 yr. They are not strictly periodic, 
since the surface density of the layered part of the disc is not 
completely smooth, as it is incessantly perturbed by the ring 
instability. At t = 700 yr the whole dead zone contained in 
the computational domain is ignited and the accretion rate 
increases up to 4.7 x 10“® M 0 yr“^. Since the mini-outburst 
is stopped by the outer boundary, and the peak accretion 
rate is probably underestimated. 


3.4 Model 4. (u = acf/Q., Oa = 10 odz = 0) 

Simulation that produced Model 1 is now repeated with the 
modified viscosity prescription. Since in the modihed pre¬ 
scription there is no direct connection between the thickness 
of the active layer and the viscosity coefficient, the ring in¬ 
stability needs much more time to develop. The growth rate 
of the instability is approximately one order of magnitude 
lower than in Model 1. By the end of the simulation only 
the innermost ring is separated from the dead zone, while 
the remaining ones are visible only as density enhancements 
within the dead zone. No ignition events are observed. The 
radial extent of the innermost ring is approximately twice 
as big as the one in Model 1, and its mass is 5 x 10“® M© 
(approximately five times more than in Model 1). 
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3.5 Model 5. {v = ac^/^l, Oa = 10 odz = 10 A 

This model differs from Model 2 in the same way as Model 4 
from Model 1, i.e. only in the formula employed to de¬ 
scribe the viscosity. With the modified formula the resid¬ 
ual viscosity in the formally dead area is high enough to 
prevent the formation rings before the dead zone ignites at 
^ 700 yr. As a result of ignition the accretion rate increases 
5 X 10“® MQyr“^ to 3 x 10~® M 0 yr“^. This mini-outburst 
involves the whole dead zone in the computational domain, 
so the peak accretion rate is probably underestimated, as 
in the case of Model 3. Approximately 5 x 10“® Mq were 
accreted onto the star during this mini-outburst. 


3.6 Model 6. [u = aCg/Q, Oa = 10 qdz = 10 

In this case, the high viscosity in the dead zone completely 
suppresses the ring instability and the surface density re¬ 
mains smooth in the layered region. Recurrent ignitions and 
replenishments of the dead zone lead to the periodic oscilla¬ 
tions of its inner boundary. The period of these oscillations 
is around 30 yr. Since the mass accumulates in the whole 
dead zone, a progressively larger part of it is ignited at each 
oscillation, and the amplitude of the oscillations generally 
increases. This increase is accompanied by a slow secular 
increase of the accretion rate at the inner boundary (from 
1 X 10“® MQyr“^ in the beginning to 1.3 x 10“® MQyr”’^ at 
364 yr). After each ignition the accretion rate temporarily 
increases by approximately 10 %. 


3.7 Model 7. (u = ac^/Q, Oa = 5 x 10 
CtDZ = 5 x 10~^) 

The model evolves similarly to Model 6 . The accretion 
rate slowly grows from 2 x 10“® in the beginning to 5.5 x 
10~® M 0 yr“^ at t = 600 yr. The relative changes during 
mini-outbursts are slightly higher, reaching approximately 
30%. At t = 400 yr a strong convection occurs in the inner 
a-disc which makes the accretion rate curve noisy. 


3.8 Model 8. {u = aCg/fl, Qa = 2 x 10 
aoz = 2 x 10 “®) 

This model exhibits very regular mini-outbursts with pe¬ 
riod ~ 150 yr. The accretion rate at the inner boundary 
changes from 2 x 10 “® in the quiet state to the peak value 
6 x 10“® M 0 yr“®. The amplitude of the oscillations re¬ 
mains constant, however, double mini-outbursts occur at 
later stages of the evolution (starting at 1300 yr). The time 
separation of the two peaks in such double mini-outburst 
is approximately 50 yr. The outer part of the dead zone is 
not affected by the oscillations - it remains stationary (see 
>14.11 . 

3.9 Summary 

In general, the numerical models exhibit two most remark¬ 
able phenomena. If aoz = 0 (models 1 and 4), the dead zone 
tends to split into rings. On the other hand, if qdz = O.loa 
(models 3 and 6 - 8 ), the dead zone remains unsplit, but the 
inner boundary of the dead zone starts to oscillate. The mass 


accumulated in this part of the dead zone is periodically de¬ 
pleted. These oscillations are accompanied by the variations 
in the accretion rate in the inner active disc. Models 2 and 
5 (where aoz = O.Olaa are the transitions between the two 
cases. 

Since the ring instability observed in models 1 and 4 
(and partially also 2 and 5) was described in detail in Pa¬ 
per I, in the following we concentrate on the description and 
explanation of the oscillations. 


4 OSCILLATIONS OF THE INNER 
BOUNDARY OF THE DEAD ZONE 


If udz a 0 some accretion in the dead zone occurs. However, 
as long as the accretion in the active layers dominates, the 
accretion rate is, in general, a function of radius. It leads to 
the accumulation of the mass and the growth of the dead 
zone surface density 2Edz at a rate 


2Sdz = 


1 dM 
2nr dr 


( 3 ) 


where M is the r-dependent accretion rate given by the mo¬ 
mentum conservation law 

M = 127rr®^^—(r'aSa-I- i^dzEdz) , (4) 

or 

where Ea, Vs. and zzdz are surface density of the active layer, 
viscosity in the active layer and in the dead zone, respec¬ 
tively. Following the standard approach of vertically averag¬ 
ing we assume the viscosities in forms 

c2 

l^a.DZ = Cla.DZ-^ (5) 

where Cs,m = is the sound speed in the mid-plane of 

the disc. 

The vertical temperature profile of a disc with z - 
dependent dissipation rate was obtained bv lHnbe~ 

In an optically thick approximation it can be written as 

T*{m) = [r(m) - re(m)] ( 6 ) 

using a mass-depth coordinate m = pdz' changing from 
0 (at the disc surface) to E/2 (at the mid-plane). r(m) is 
the optical depth and 


rg{m) = / K{m')9{m')dm' 


is the viscosity-weighted optical depth, where 


( 7 ) 


9{m) = 


v(m')dm' 
iZaEa -I- ZZdzEdZ 


( 8 ) 


is the viscosity weight - a monotonically increasing function 
between 9(0) = 0 and 9(E/2) = 1. In our model it assumes 
the form 



_ tzam _ 

UaSa + i^DZ^DZ 
t^aSa + t^Dz(77i-Sa) 
i^aSa + i^DZ^DZ 


for ( 0 , Ea) 
for (Ea, E/2) 


( 9 ) 


The effective temperature, Te, is given by the formula 
Ts = -7—^^ (EaZZa -|- EoZi^Dz) , (10) 
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Model 1 


Model 4 


Model 5 


Model 2 


Model 3 


Models 


Figure 1. Models 1-6: the top figures show the evolution of the disc surface density (levels of grey), the dashed line denotes the boundary 
between the dead zone and the active parts of the disc in the mid-plane. The bottom figures show the accretion rate at the inner boundary 
of the computational domain. 
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Figure 2. Models 7 and 8. As in Fig. 1. 


where a is the Stefan-Boltzmann constant, is the Kep- 
lerian angular velocity. Using Eq. 0-® and 0 we can 
express the mid-plane temperature Tm = r(E/2) as 


„4 _ 3 4 QaSa -|- aDzEDz(EDZ + 2Ea) 

“ ^ 8 " aaEa anzEnz 


( 11 ) 


Further, we assume that the opacity is a piecewise 
power-law function of temperature, k = kqT^, where the 
coefficients Kn and b f or dif ferent temperature ranges are 
taken from iBell fc LirJ il994l . Inserting Eq. itm into 0 we 
obtain a formula for the accretion rate of the layered disc 



A|nT^ri/"(aaEa-faDzEDz) (12) 

[aaEa-I-anzEnz (Edz + 2Ea)] I , 


where r = kE/ 2 is the optical depth at the mid-plane. In the 
same solution the effective temperature is explicitly given by 

Te,c. = . (15) 

Combining Eqs. iHH and ilT^ with the opacity formula we 
get 


/ 1280CT^mH \ 5/4 

V 27kBa,n J ■ 


(16) 


The inner boundary of the dead zone oscillates for rea¬ 
sons explained Fig. El Eign and Emtn intersect at Ri dehned 
by Gammie as the inner boundary of the dead zone - radius 
at which limit layered disc with Edz = 0 has the mid-plane 
temperature Tum 


Ri = 


27kBQ.a. 

320cr/imH 


2/3 


(17) 


which will be used in am 

As more and more mass accumulates in the dead zone, 
the mid-plane temperature of the disc grows and when it 
reaches Tum the disc is made active. We can estimate the 
surface density at the moment of ignition, Eign, by substi¬ 
tuting Tiim for Tm in Eq. Illll and using a metal-grains dom¬ 
inated opacity with kq = 0.1 and b — 0.5 (appropriate for 
temperatures between 203 K and 2290p^^^® K). We get 


/ 320a- /rrwH 5/2 
V 27 kB^aoz 


Qa — Qpz ^2 
ClDZ 


, 1/2 


(13) 


On the other hand, there is a minimum surface density 
Emtn necessary to maintain the disc active. It is the surface 
density at which the mid-plane temperature of a fully active 
disc (i.e. of a disc in which the dead zone vanishes) is equal 
to Tiim. We estimate it with the help of the standard a-disc 
solution, in which the mid-plane temperature Tm,c« and the 
effective temperature Te,a obey the relation 



(14) 


In numerical models, however, the inner boundary of the 
dead zone is shifted outwards with respect to Ri. This is 
because of two effects that increase the temperature in the 
inner part of the layered disc: the radiation tran sfer from 
the a djacent hotter active disc (see also footnote in lGamrnI3 
1199611 : and the residual viscous dissipation in the dead zone. 

As soon as the surface density E exceeds Eign, the dead 
zone is ignited. The ignition always occurs at the inner 
boundary of the dead zone, because Edz is the highest there 
(it can be shown by solving Eq.Elnumerically) and Eign is an 
increasing function of the radius. Therefore, the dead zone 
is ignited from inwards. 

In the ignited area the increasing viscosity causes tem¬ 
perature and accretion rate to grow. The radiation flux from 
the ignited region heats up and ignites the adjacent annu¬ 
lus of the dead zone, so that the active area spreads out¬ 
wards. On the other hand, the increasing accretion rate 
tends to decrease the surface density of the ignited region, 
since the mass is transported from the outer disc at a con¬ 
stant rate independent on whether the dead zone has ignited 
or not. When E drops below Emtn the ignited region becomes 
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Figure 3. The mechanism of the oscillations of the inner bound¬ 
ary of the dead zone. Top: the mid-plane temperature of the nu¬ 
merical model compared to the analytical one calculated for a disc 
with the empty dead zone (i.e. Edz = 0)- Bottom: the profiles of 
Sign a^nd Smtn- During mini-outbursts the surface density in the 
inner part of the dead zone oscillates between these two curves. 



r[AU] 

Figure 5. The stationary profiles of the surface density for dif¬ 
ferent accretion rates at the inner boundary. The sharp changes 
of slopes are due to switching between different opacity regimes 
at temperatures 203 K and 167 K. The dashed lines show Sign 
and Smtn, the horizontal dotted line denotes a value of 2Sa. The 
distinctive points Ri, Ri and R 2 for M = 10”^ Mq yr“^ model 
are denoted by the vertical dotted lines. 



Figure 4. A time-sequence of surface density profiles during the 
mini-outburst in model 6 that occurs between 28 and 53 yr. Solid 
lines (28 — 38.7 yr): the boundary of the dead zone receding from 
the star; dot-dashed line (42.8 yr) the boundary has reached the 
maximum distance from the star; dashed lines (46.4 — 52.7 yr): 
the boundary approaches the star. 

dead, the layered structure is reestablished, and the accu¬ 
mulation of mass starts again. This process is illustrated in 
Fig. a which shows a time sequence of surface-density pro¬ 
files during one such ’’mini-outburst” observed in the nu¬ 
merical model. As one might expect, this behaviour is qual¬ 
itatively similar to changes in the surface density profile of 
disc s in cataclysmic variables (CVsi during an outburst; see 
e.g. iMineshige fc Osakil (^8^. However, while in CVs two 
types of outbursts occur (outside-in and inside-out), only 
the latter seems to be possible in layered discs. As we al¬ 
ready mentioned, this is because at smaller radii the surface 
density of the dead zone increases faster. In principle, it is 
conceivable that an outburst might produce a density pro¬ 
file with which the next ignition would occur away from the 
inner edge. We cannot entirely exclude such possibility, al¬ 
though in our simulations, covering a few hundred outbursts 
in total, such a case has never been observed. 


4.1 Stationary solution 

The surface density of the dead zone evolves according to 
Eq. 0. As Edz grows and the accretion in the dead zone 
becomes more and more important, the accumulation rate 
Edz decreases. At larger radii a stationary state with radi¬ 
ally constant accretion rate may be reached. 

We search for the stationary surface density profile Sstat 
by solving Eq. 0 numerically with a 1- D co de similar 
to tha t described in lArmitage et al.l i2nnill and ISteninsO 
lIlQQlil . The computational domain (0.1 — 10 AU) is divided 
into 500 cells equidistantly, and the surface density E is set 
to some initial value (e.g. 2Ea). Then, the the mid-plane 
temperature is calculated from Eq. o for every grid cell 
and it is checked if the cell belongs to the layered part of 
the disc (if Tm < Tiim and E > 2Ea) or to the standard 
a-disc (otherwise). In the the later case, the temperature is 
recomputed according to Eq. itni to obtain the consistent 
a-disc solution there. 

The opacity coefficients kq and b are determined using 
the first three Bell & Lin regimes only since the analytical 
solution does not always exist for the remaining ones (i.e. we 
take into account ice dominated opacity with kq = 2 x 10”“^ 
and 6 = 2 below 167 K, the jump due to ice evaporation 
with Ko = 2 X 10^® and 6 = —7 between 167 and 203 K, and 
metal grains dominated opacity with kq = 0.1 and 6 = 0.5 
between 203 and 2290^^^"*® K). Our opacity formula may 
be incorrect for the part of the inner a-disc (and, occasion¬ 
ally, also for the several cells of the layered region where the 
temperature exceeds 2290p^^’*® K). However, using this toy 
model we are not particularly interested in the very inner 
disc, which is much better described by the TRAMP sim¬ 
ulation anyway. Moreover, recent works on opa city in pro¬ 
toplanetary discs (see e.g. ISemenov et al ]|20 q 1) show that 
the metal grains evap orate at a s l ightly higher temperature 
than that adopted bv iBell fc LirJ lll99.^ . 

Subsequently, the accretion rate M is computed at the 
boundaries of each cell from Eq. JUJ for the layered disc or 
from 
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Ma = (18) 

for the a-disc. Finally, the amount of mass proportional to 
some small time-step dt is transported between cells and E 
is updated. The accretion rate Mob at the outer boundary 
of the domain is a free parameter, and a free outflow con¬ 
dition is imposed at the inner boundary. The computation 
ends when the accretion rate throughout the computational 
domain converges to Mob, i.e. when the stationary state 
is reached. Typically, it happens after several times 10® yr 
depending on Mob- 

The stationary layered disc surface density profiles Estat 
for M = 10“®, 10“’’ and 10“® are shown in Fig. They 
end at the radius Ri where they cross Eign curve, since the 
stationary layered disc solution does not exist below Ri. 
The radius R 2 at which Estat curve crosses Emtn is another 
distinctive point: it separates the region of the dead zone 
which after the ignition is able to maintain the temperature 
above Tiim from the one which is not. The dead zone ends 
at radius Ri where the surface density drops below 2Ea. 

In this manner, the dead zone can be divided into three 
parts: a) the non-stationary periodically depleted and re¬ 
plenished part between Ri and Ri, b) the stationary, but 
"combustible” part between Ri and R 2 , and c) the station¬ 
ary "incombustible” part between R 2 and Ro- 

An upper limit for the mass accreted during one mini¬ 
outburst in region a) can be obtained by integration the 
difference between Eign and Emtn curves there 

/•Rl 

Ma = 27v / r(Sign - Tlnitn)dr . (19) 

J Ri 

Taking values Ri from the numerical model R\ = 0.15, 0.38 
and 0.96 AU for M = 10“®, 10“’’ and 10“® Mq yr“’, respec¬ 
tively, we get 

{ 7.2 X 1O“'’M0 for M = 10“® Mq yr“’ 

2.6 X 1O“®M0 for M = 10“'^ M 0 yr“’ (20) 

5.2 X 1O“'‘M0 for M = 10“® M 0 yr“’ 

The amount of mass available for the accretion in the 
region b) may be interesting in relation to FU Orionis type 
outbursts. It can be determined as 
rR2 

Mh = 2'K / r(Estat - Emtn)dr . (21) 

J Ri 

Integrating numerically the stationary solutions obtained by 
1-D model we get 

{ 2.0 X 10“®Mq for M = 10“® Mq yr“’ 

2.7 X 10“^Mq for M = 10“’’ Mq yr“’ (22) 

O.OIMq for M = 10“® Mq yr“’ 

The typical amount of mas s accreted during the FU Or i 
type outburst is ~ 0.01 Mq dllartmann fc KenvnullToOfil 
which is comparable to the mass stored and available for 
the accretion in the dead zone in the case with M = 
10“® M0yr“’. 


5 DISCUSSION 

The two-dimensional radiation-hydrodynamic simulations of 
the layered disc we carried out showed that the behaviour of 
the disc strongly depends on the presence of a small residual 


viscosity in the dead zone q:dz. In the models with qdz = 0 
the ring instability develops as described in Paper I. The 
rings may eventually decay which would temporarily in¬ 
crease the accretion rate onto the central star by a factor 
of several. On the other hand, in the models with odz set to 
10% of the viscosity in the active parts of the disc, the dead 
zone remains contiguous, but its inner boundary oscillates. 

These oscillations are more or less regular with period 
between several tens and several hundreds years depending 
on model parameters. They can be accompanied by varia¬ 
tions of the accretion rate on a time-scale as short as ten 
years, at a relative amplitude reaching from a few per cent 
to 200 - 300 per cent in the high viscosity Model 8. These 
variations would be rather difficult to observe, however it 
does not seem absolutely impossible. Since the properties of 
these oscillations are quite sensitive to physical conditions 
in the disc, such observations could help to constrain the set 
of models of protoplanetary discs. 

Similar variati ons of the accretion rate were also found 
bv ISteoinskil ^1999^ using an one-dimensional hydrodynamic 
code. However, the period was much longer than in our 
model (10® — lO"’ yr). The reason for that may be that in 
his models the ignition of the inner part of the dead zone 
is much more difficult due to the absence of the radiation 
transport (and, generally speaking, due to oversimplified de¬ 
scription of the interactions between the inner hot a-disc and 
the layered region). 

We give an analytical description of these oscillations 
based on the vertically averaged model modified for the case 
of layered disc. It illustrates the mechanism of the dead zone 
oscillations and tests the numerical model. 

Finally, we perform a simple 1-D model of the layered 
disc based on the analytical solution to follow the evolution 
of the whole layered region for a long time. We found that 
with Qdz yf 0 a stationary layered region may exist at higher 
radii (this was also observed in numerical Model 8). Part 
of this stationary region is still combustible, i.e. it may be 
made active by increasing the temperature by some external 
event, and part of its mass is available for accretion. 

We compare the amount of mass stored in the com¬ 
bustible region to the mass accreted during a FU Ori type 
outburst following the idea that the dead zone can serve as a 
mass reservoir for such event. We found, that these masses 
are comparable for the case of massive disc with the high 
accretion rate, which is in agreement with the observation 
that the FU Ori objects are very young systems. 
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